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SOME  SCHWARZ  METHODS  FOR  SYMMETRIC  AND 
NONSYMMETRIC  ELLIPTIC  PROBLEMS 

OLOF  B.  WIDLUND  * 

Abstract.  This  paper  begins  with  an  introduction  to  additive  and  multipHcative  Schwarz  methods. 
A  two-level  method  is  then  reviewed  and  a  new  result  on  its  rate  of  convergence  is  established  for  the 
case  when  the  overlap  is  small.  Recent  results  by  Xuejun  Zhang,  on  multi-level  Schwarz  methods, 
are  formulated  and  discussed.  The  paper  is  concluded  with  a  discussion  of  recent  joint  results  with 
Xiao-Chuan  Cai  on  nonsymmetric  and  indefinite  problems. 

Key  Words,  domain  decomposition,  Schwarz  methods,  finite  elements,  nonsymmetric  and  indefi- 
nite elliptic  problems 
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1.  Introduction.  Over  the  last  few  years,  a  general  theory  has  been  developed 
for  the  study  of  additive  and  multiplicative  Schwarz  methods.  Many  domain  decom- 
position and  certain  multigrid  methods  can  now  be  successfully  analyzed  inside  this 
framework.  Early  work  by  P.-L.  Lions  [23],  [24]  gave  an  important  impetus  to  this 
effort.  The  additive  Schwarz  methods  were  then  developed  by  Dryja  and  Widlund 
[15],  [16],  [17],  Matsokin  and  Nepomnyaschikh  [25]  and  Nepomnyaschikh  [26]  and  oth- 
ers. Recent  efforts  by  Bramble,  Pasciak,  Wang  and  Xu  [4]  and  Xu  [36]  have  extended 
the  general  framework  making  a  systematic  study  of  multiplicative  Schwarz  methods 
possible.  The  multiplicative  algorithms  are  direct  generalizations  of  the  original  al- 
ternating method  discovered  more  than  120  years  ago  by  H.  A.  Schwarz  [30].  We  note 
that  most  of  the  work  in  recent  years  has  focused  on  the  positive  definite,  symmetric 
case. 

While  this  theory  is  quite  general,  the  applications  so  far  have  primarily  been  to 
the  solution  of  the  often  large  linear  systems  of  algebraic  equations,  which  arise  in 
the  finite  element  discretization  of  elliptic  and  parabolic  boundary  value  problems. 
As  shown  in  P.-L.  Lions  [24],  the  classical  Schwarz  algorithms  can  conveniently  be 
described  in  terms  of  subspaces  of  the  given  space.  The  relevant  error  propagation 
operator  of  a  particular  Schwarz  method  can  be  written  as  a  polynomial  of  orthog- 
onal projections  onto  these  subspaces.  The  use  of  these  projections  in  computations 
involves  the  evaluation  of  the  residual  of  the  original  finite  element  problem  and  the 
exact,  or  approximate,  solution  of  several  finite  element  problems  on  subregions.  An 
additional  coarse  discrete  model  is  also  often  used  to  enhance  the  rate  of  convergence. 
For  a  discussion  of  many  applications,  see  Dryja  and  Widlund  [17].  For  other  current 
projects,  which  also  use  the  Schwarz  framework,  see  Dryja  and  Widlund  [18],  [20]  and 
Dryja,  Smith  and  Widlund  [14]. 


'  Courant  Institute  of  Mathematical  Science,  New  York  University,  251  Mercer  Street,  New  York, 
N.Y.  10012.  This  work  was  supported  in  part  by  the  National  Science  Foundation  under  Grant  NSF- 
CCR-8903003  and  in  part  by  the  U.S.  Department  of  Energy  under  contract  DE-FG02-88ER250553 
at  the  Courant  Mathematics  and  Computing  Laboratory. 
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This  paper  is  organized  as  follows.  In  Section  2,  we  first  briefly  review  the  theory 
for  the  positive  definite,  symmetric  case  and  then,  to  provide  an  example,  turn  to 
a  special  case  already  discussed  in  our  first  papers  on  Schwarz  methods,  cf.  Dryja 
and  Widlund  [15],  [16].  This  additive  method  uses  two  levels  of  discretization  corre- 
sponding to  a  coarse  space  and  a  number  of  local  finite  element  spaces.  The  former 
plays  a  role  similar  to  the  coarse  mesh  solver  in  a  multigrid  method;  it  provides  global 
transportation  of  information  across  the  region. 

In  our  original  proof  that  the  condition  number  for  this  Schwarz  method  is  inde- 
pendent of  the  number  of  subspaces  and  the  mesh  size,  it  is  important  to  assume  that 
there  is  a  relatively  generous  overlap  between  the  subregions.  However,  numerical 
experiments  have  shown  that  such  methods  often  perform  quite  satisfactorily  with 
Uttle  overlap,  cf.  Bj0rstad,  Moe  and  Skogen  [1],  Bj0rstad  and  Skogen  [2],  Cai  [7], 
[8]  and  Cai,  Gropp  and  Keyes  [9].  In  Section  3,  we  prove  a  theorem  on  additive 
Schwarz  methods  with  minimal  overlap  to  help  provide  an  explanation.  This  result 
has  recently  been  obtained  in  joint  work  with  Maksymilian  Dryja. 

In  Section  4,  we  describe  a  recent  result  by  Xuejun  Zhang  [37],  a  Courant  Institute 
student,  on  multilevel  additive  Schwarz  methods.  He  has  shown  that  the  condition 
number  of  a  family  of  such  methods  grows  at  most  linearly  with  the  number  of  levels. 
The  parallel  multilevel  preconditioner  of  Bramble,  Pasciak  and  Xu  [5]  belongs  to  this 
family  and,  consequently,  the  estimate  of  the  condition  number  of  the  BPX  operator 
can  be  improved.  For  earlier  results  and  a  general  discussion  of  multilevel  methods, 
see  Dryja  and  Widlund  [19]  and  Xu  [35]. 

In  Section  5,  we  discuss  recent  joint  work  with  Xiao-Chuan  Cai  on  nonsymmet- 
ric  and  indefinite  elliptic  problems.  Cai's  dissertation  [7]  was  a  pioneering  study  on 
nonsymmetric  problems.  He  discovered  that  the  use  a  coarse  model  is  even  more 
important  than  in  the  positive  definite,  symmetric  case.  Thus,  already  in  the  non- 
symmetric,  positive  definite  case,  the  use  of  a  sufficiently  fine  coarse  mesh  model  is 
required  in  order  for  the  additive  Schwarz  operator  to  have  a  positive  definite,  sym- 
metric part.  Confining  the  spectrum  of  the  additive  Schwarz  operator  to  the  right 
half  plane  allows  us  to  derive  strong  results  on  the  convergence  rate  of  GMRES,  a 
conjugate  gradient-like  method,  which  is  effective  for  many  nonsymmetric  problems, 
cf.  Saad  and  Schultz  [28].  Generally,  it  is  known  that  the  performance  of  iterative 
methods  of  the  conjugate  gradient  family  suffers  if  the  origin  of  the  complex  plane 
is  surrounded  by  eigenvalues.  The  experimental  evidence  is  also  very  clear;  the  per- 
formance is  clearly  enhanced  in  the  presence  of  a  coarse  space  with  sufficiently  many 
degrees  of  freedom,  cf.  Cai  [7],  Cai,  Gropp  and  Keyes  [9]  and  Cai  and  Widlund  [11]. 

Our  joint  work  concerns  both  additive  and  multiplicative  algorithms,  cf.  Cai 
and  Widlund  [11],  [10].  In  the  second  paper,  the  theory,  previously  developed  for 
multiplicative  methods  by  Bramble  et  al.  [4]  for  the  positive  definite,  symmetric  case, 
is  generalized  to  a  class  of  nonsymmetric  and  indefinite  elliptic  problems. 
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2.   An  introduction  to  Schwarz  methods.  As  usual,  we  write  our  continuous 
and  finite  element  elliptic  problems  as 

a{u,v)  =  f{v),  V  t;  G    F  , 

and 

(1)  a{u,,VH)  =  f{vH),yv,e    V\ 

respectively.  We  assume,  in  the  next  few  sections,  that  the  problem  is  selfadjoint  and 
elliptic  in  a  suitably  chosen  space  V  and  that  a{u,v)  is  bounded  in  V"  x  V.  In  the  case 
of  Poisson's  equation,  the  bilinear  form  is  defined  by 


Vu  •  Vu  dx 


To  avoid  unnecessary  complications,  we  first  confine  our  discussion  to  this  operator,  to 
homogeneous  Dirichlet  conditions  and  to  continuous,  piecewise  linear  finite  elements. 

In  the  classical  formulation  of  Schwarz's  alternating  method,  two  overlapping 
subregions,  Q[  and  Q'2  are  used;  the  union  of  the  two  is  fi.  There  are  two  sequential, 
fractional  steps  in  which  the  approximate  solution  of  the  elliptic  equation  on  H  is 
updated  by  solving  the  problem  restricted  to  one  of  the  subregions,  one  at  a  time. 
The  most  recent  values  of  the  solution  are  used  as  boundary  values  on  the  part  of 
dQ,[,  the  boundary  of  Q(,  which  is  not  a  part  of  dQ. 

The  finite  element  version  of  the  algorithm  can  conveniently  be  described  in  terms 
of  projections  P,  :  V''  —y  F,'',  defined  by 

(2)  a(P,Vh,  M  =  a{vh,  <f>h),  ^4>h  e  V,''  . 

Here  V'^  is  the  finite  element  space  on  Q,  and  V/^  =  HQ{Q,'-)f]V'^.  It  is  easy  to  show 
that  the  error  propagation  operator  of  this  multiphcative  Schwarz  method  is 

(I  -  P,){I  -  P,). 
This  algorithm  can  therefore  be  viewed  as  a  simple  iterative  method  for  solving 

{Pi  +  P2  -  P2Pl)Uh  =  9h, 

with  an  appropriate  right-hand  side  gh.- 

This  operator  is  a  polynomial  of  degree  two  and  thus  not  ideal  for  parallel  com- 
puting since  two  sequential  steps  are  involved.  This  effect  is  further  pronounced  if 
more  than  two  subspaces  are  used  even  if  the  degree  of  the  polynomial  representing 
the  multiplicative  algorithm  often  is  lower  than  maximal  because  the  product  of  two 
projections  associated  with  subregions,  which  do  not  overlap,  vanishes. 

We  note  that  it  is  often  advantageous  to  group  subspaces,  which  do  not  intersect, 
into  sets.  The  subspaces  of  each  set  can  then  be  merged  and  the  value  of  A^  reduced. 
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The  decrease  in  the  number  of  fractional  steps  makes  the  algorithm  easier  to  paral- 
lelize. Numerical  experiments  also  show  that  the  convergence  rate  often  is  enhanced 
if  the  number  of  subspaces,  created  in  this  way,  is  small. 

In  the  additive  form  of  the  algorithm,  we  work  with  the  simplest  possible  poly- 
nomial of  the  projections:  The  equation 

(3)  PuH  =  {Pi  +  P2  +  ---  +  PN)uk  =  g'k  , 

is  solved  by  an  iterative  method.  The  choice  of  this  polynomial  can  also  be  motivated 
by  first  introducing  a  relaxation  parameter  r  in  the  multiplicative  algorithm,  which, 
for  the  case  oi  N  =  2,  results  in  the  polynomial  r(Pi  +  P2  —  tP^Pi).  CanceUing  the 
factor  r  and  letting  r  decrease  to  zero,  we  again  obtain  the  first  order  polynomial  of 
equation  (3). 

Since  the  operator  P  is  positive  definite,  symmetric,  with  respect  to  the  bilinear 
form,  the  iterative  method  of  choice  is  the  conjugate  gradient  method.  Equation  (3) 
must  have  the  same  solution  as  equation  (1),  i.e.  the  correct  right-hand  side  must 
be  found.  Since  by  equation  (1),  a{uh,<t>h)  =  f{<f>h)  ,  the  right-hand  side  g'f^  can  be 
constructed  by  solving  equation  (2)  for  all  values  of  i  and  adding  the  results.  Similarly 
the  operator  P  of  equation  (3)  can  be  applied  to  any  given  element  of  V'^  by  applying 
each  projection  P,  to  the  element.  Most  of  the  work,  in  particular  that  which  involves 
the  individual  projections,  can  be  carried  out  in  parallel. 

We  now  describe  the  special  additive  Schwarz  method  introduced  in  Dryja  and 
Widlund  [15],  cf.  also  [16]  and  Dryja  [13].  We  start  by  introducing  two  triangulations 
of  Q.  into  nonoverlapping  substructures  Q,  and  elements.  We  obtain  the  elements  by 
subdividing  the  substructures.  We  always  assume  that  these  triangulations  are  shape 
regular,  cf.  e.g.  Ciarlet  [12].  In  this  additive  method,  we  use  overlapping  subregions 
obtained  by  extending  each  substructure  to  a  larger  region  Q'-.  We  first  assume  that 
the  overlap  is  generous  with  the  distance  6,  between  the  boundaries  dQ,  and  5f2J 
bounded  from  below  by  a  fixed  fraction  of  H,,  the  diameter  of  Q,.  We  also  assume 
that  d^\  does  not  cut  through  any  element.  We  make  the  same  construction  for  the 
substructures  that  meet  the  boundary  except  that  we  cut  off"  the  part  of  f)'  that  is 
outside  of  Q. 

Our  finite  element  space  is  represented  as  the  sum  of  N-(-l  subspaces 

V'  =  v^'  +  vl'  + ...  +  vl^ . 

The  first  subspace  Vq  is  equal  to  V" ,  the  space  of  continuous,  piecewise  linear  func- 
tions on  the  coarse  mesh  defined  by  the  substructures  Q.i.  The  other  subspaces  are 
related  to  the  subdomains,  in  the  same  way  as  in  a  traditional  Schwarz  algorithm,  i.e. 

It  is  often  more  economical  to  use  approximate  rather  than  exact  solvers  of  the 
problems  on  the  subspaces.  We  use  the  following  notations.  bi{u,v)  is  an  inner 
product  defined  on  V^''  x  V;''.  We  assume  that  there  exists  a  constant  uj  such  that 

(4)  a{u,u)  <ujb,{u,u)  yueV,^  . 
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An  operator  T,  :  V''  -^  V/",  which  replaces  P.,  is  defined  by 

(5)  b,{T,u,  </..)  =  a{u,  <f>h),  y<f>h  6  V,''  . 

It  is  easy  to  show  that  the  operator  T,  is  positive  semidefinite  and  symmetric  with 
respect  to  a(-,  •)  and  that  the  minimal  constant  u>  in  equation  (4)  is  ||r,||a. 

In  order  to  estimate  the  rate  of  convergence  of  our  special  additive  Schwarz 
method,  or  any  other,  we  need  upper  and  lower  bounds  for  the  spectrum  of  the 
operator  relevant  in  the  conjugate  gradient  iteration.  The  technique  used  to  obtain 
lower  bounds  is  often  associated  with  P.-L.  Lions,  cf.  [24].  One  form  of  this  result  is 
given  in  the  following  lemma. 

Lemma  1.  Let  T,  be  the  operator  defined  in  equation  (5)  and  let  T  =  To  +  Ti  + 
\-Tn.  Then 

a{T~^u,u)  =    min   ^6,(u,, u,),  u,  G  V,  . 

Therefore,  if  a  representation,  u  =  X^u,,  can  be  found,  such  that 

^6,(u,,u,)  <  Coa(u,ti),    WeV^, 
then 

The  upper  bound  for  the  spectrum  of  T  can  often  conveniently  be  obtained  in  terms 
of  strengthened  Cauchy- Schwarz  inequalities  between  the  different  subspaces.  Note 
that  we  now  exclude  the  index  0;  the  coarse  subspace  is  treated  separately. 

Definition  1.  The  matrix  £  =  {£,_,}  is  the  matrix  of  strengthened  Cauchy- 
Schwarz  constants,  i.e.  Sij  is  the  smallest  constant  for  which 

(6)  \a{vi,Vj)\  <  e,j\\vi\\a\\vj\\a  ,Vvi  €  V,  ,Vdj  €  Vj  ,ij  >  1  , 

holds. 

The  following  lemma  is  easy  to  prove. 

Lemma  2.  Let  p{£)  be  the  spectral  radius  of  the  nnatrix  £.  The  operator  T  satisfies 
the  following  upper  bound, 

T<uj{p(£)  +  l)I  . 

Combining  these  two  lemmas,  we  obtain 

Theorem  1.  The  condition  number  of  the  operator  of  the  additive  Schwarz 
method  satisfies 

K{T)<u{p{£)  +  \)Cl. 

In  the  special  case  considered  in  this  section,  it  is  easy  to  show  that  there  is  a 
uniform  upper  bound;  see  Section  3,  where  an  estimate  of  Co  is  also  given. 
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It  is  clear  that  Co  does  not  increase  if  we  expand  (or  add)  individual  subspaces. 
This  follows  from  the  fact  that  there  is  a  larger  choice  in  selecting  u,  G  K-  If  we  can 
expand  the  subspaces  without  worsing  the  upper  bound,  and  that  is  often  possible, 
our  estimate  of  the  condition  number  k{T)  can  improve.  On  the  other  hand,  a  larger 
subspace  also  means  that  the  subproblems  have  more  variables  and  that  they  are 
worse  conditioned. 

In  the  multiplicative  case,  we  need  to  provide  an  upper  bound  for  the  spectral 
radius  of  the  error  propagation  operator 

(7)  Ej  =  {I-Tj)---{I-To). 

The  following  theorem  is  a  variant  of  the  important  results  of  Bramble,  Pasciak,  Wang 
and  Xu  [4]  and  Xu  [36]. 

Theorem  2.  In  the  symmetric,  positive  definite  case, 


p{Ej)  < 


(2-c^) 


\        {2u^p{£r  +  l)CS 


We  note  that  this  formula  is  useless  if  a;  >  2.  Since  ||/  —  T,||a  >  1  if  \\T,\\a  >  2,  the 
assumption  that  u;  <  2  is  most  natural. 

3.  The  case  of  a  small  overlap.  As  we  have  already  pointed  out,  the  arith- 
metic work  grows  with  the  overlap  while,  at  the  same  time,  the  smallest  eigenvalue  of 
the  operator  P,  or  T,  can  increase,  often  leading  to  a  decrease  of  the  condition  num- 
ber. Numerical  experiments  suggest  that  it  is  often  advantageous  to  use  a  minimal 
overlap,  with  a  single  layer  of  interior  mesh  points  common  to  neighboring  subregions, 
if  we  wish  to  minimize  the  running  time  of  the  algorithm  rather  than  the  number  of 
iterations,  cf.  Bj0rstad,  Moe  and  Skogen  [1],  Bj0rstad  and  Skogen  [2],  Cai  [7],  [8]  and 
Cai,  Gropp  and  Keyes  [9].  Some  experiments  also  indicate  that  for  large  values  of 
H/h,  the  condition  number  decreases  by  about  a  factor  of  two  when  the  overlap  is 
doubled.  All  these  results  are  for  problems  in  the  plane.  For  a  discussion  of  a  similar 
phenomenon  with  a  related  algorithm,  see  also  Smith  [33]  or  Chapter  4  of  Smith  [32]. 
It  is  important  to  note  that  in  these  experiments,  and  in  the  proof  given  below,  the 
special  coarse,  global  subspace  is  very  important.  Without  it,  these  algorithms  can 
be  very  slow. 

We  will  give  an  estimate  which  helps  explain  the  results  of  these  experiments. 
Our  proof  is  valid  in  two  as  well  as  three  dimensions  and  this  suggests  that  additive 
Schwarz  methods  with  small  overlap  might  produce  satisfactory  results  even  in  three 
dimensions,  in  particular,  if  the  value  of  H/h,  and  thus  the  number  of  degrees  of 
freedom  for  each  subproblem,  is  not  enormous.  We  also  note  that  parallel  computing 
practitioners  have  expressed  a  clear  preference  for  domain  decomposition  methods 
where  the  subregions  only  share  one  layer  of  points,  cf.  e.g.  Smith  [31],  [34].  Iterative 
substructuring  methods  have  that  property  but  so  do  the  additive  Schwarz  methods 
with  minimal  overlap.   From  now  on,  C  will  denote  a  generic  constant.   To  simplify 
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the  proof,  we  assume  that  the  maximum  width  of  F^,,,  the  neighborhood  inside  the 
boundary  5Q,  that  is  in  common  with  other  subregions,  is  bounded  by  C(S,. 

Theorem  3.  In  the  case  of  exact  solvers  for  the  suhyrohlems,  the  condition 
number  of  the  additive  Schwarz  method  satisfies 

k(P)  <C(1  +  Cmax^). 

Here  Hi  is  the  diameter  of  the  substructure  Q,  and  6,  the  distance  between  the  bound- 
aries of  Q,  and  0,'-.  The  constants  are  independent  of  the  parameters  Hi,  h  and  Si. 

Proof.  The  proof  is  a  refinement  of  a  result  first  given  in  Dryja  and  Widlund 
[15];  cf.  [16]  for  a  better  discussion.  We  first  note  that  a  constant  upper  bound  for 
the  spectrum  of  P  can  be  obtained  as  follows.  For  «  >  1,  P,  is  also  an  orthogonal 
projection  of  H^{Q,[)  onto  VJ.  Therefore, 

a(P,Uh,Uk)  <  aQ'(uh,Uh). 

Since  there  is  an  upper  bound,  m,  on  the  number  of  subregions  to  which  any  x  G  f2 
can  belong,  we  have 

N 

!  =  1 

In  addition,  we  use  the  fact  that  the  norm  of  Pq  is  equal  to  one  and  obtain 

Amax(P)<(m  +  l). 

The  lower  bound  is  obtained  by  using  Lemma  1.  A  natural  choice  of  uq  is  the 
L2-projection  Qhuh  of  Uk  onto  V^ .  It  is  well  known,  cf.  e.g.  Bramble  and  Xu  [6], 
that  Qh  is  bounded  not  only  in  L2  but  also  in  H^  and  that  there  exists  a  constant, 
independent  of  h  and  i7,  such  that 

(8)  \\QHUh-Uh\\L2    <C  H\\uh\\a- 

Let  Wh  =  u/i  —  QnUh  and  let  u,  =  Ih{9,Wh)  , «  —  I,-  ■  ■  ,N.  Here  Ik  is  the  interpolation 
operator  onto  the  space  V''  and  the  9,  define  a  partition  of  unity.  These  functions  can 
be  chosen  in  many  ways,  in  particular,  as  nonnegative  elements  of  V''.  Then,  all  the 
values  of  6i  are  determined  by  those  at  the  mesh  points.  In  the  interior  part  of  Q(-, 
which  does  not  belong  to  any  other  subregion,  6i  =  1.  This  function  must  also  drop 
to  0  over  a  distance  on  the  order  of  ^,.  We  denote  the  region,  where  there  is  overlap, 
by  Ts  and  we  also  use  the  notation  F^,,  =  F^  fl  Q,.  It  is  easy  to  construct  a  partition 
of  unity  with  0  <  ^,  <  1  and 

iv.,i  <  £ . 

In  order  to  use  Lemma  1,  we  now  estimate  a(u,, u,)  in  terms  of  a(wh,Wh).  We 
consider  the  contribution  from  one  substructure  at  a  time  and  note  that,  trivially, 

an,\r«.,("M".)  =  aQ,\rs,,{'^h,Wh)  • 
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Let  K  be  an  element  in  F^,,.  Then,  using  the  definition  of  u,, 

where  6i  is  the  average  of  6i  over  the  element  K.  Using  the  fact  that  the  diameter  of 
K  is  on  the  order  of  h  and  the  bound  on  V^,,  we  obtain 

^  C  ..  ,,2 

aK{ui,  u,)  <  2ah-{wh,  Wk)  +  -p\\wh\\mK)  ■ 

To  complete  the  proof,  we  need  to  estimate  ||u^/i||£,2(r    )•  This  estimate  is  based 
on  the  following  lemma. 
Lemma  3. 

\H\Ur„)  <  CSfiil  +  |^)an.(u, u)  +  -^||u||i.(o.))  . 

Proof.    All  the  essential  ideas  of  the  proof  can  be  illustrated  by  considering  a 
square  region  (0,  H)  x  (0,  H).  Then,  since 


fy  du(x, t) 

u{x,0)  =  u{x,y)  ~ dr  , 

Jo        oy 

we  find,  after  some  manipulations,  that 

H         \u(x,0)\^dx  <2         /     \u{x,y)\^dxdy  +  H^         /     l^l^dxdy  . 
Jo  Jo    Jo  Jo    Jo      oy 

Therefore, 

Now  consider  the  integral  over  a  narrow  subregion  next  to  one  of  the  sides  of  the 
square.  Using  similar  arguments,  we  obtain 

/     [\u{x,y)\^dxdy<6^aQ,{u,u)  +  26f     \u{x,0)\'^dx  . 
Jo    Jo  Jo 

By  combining  this  with  the  previous  inequality,  we  obtain 

rH      rS  0 

Jo    Jo   l"(^'y)l^^^^2/ <  <^^an.(",")  +  2<!»(-^||u||^2(fj_j  +  /ran.(u,u))  , 


as  required. 

The  modifications  necessary  to  cover  the  case  of  an  arbitrary,  shape  regular  sub- 
structure are  routine.  | 

We  now  apply  Lemma  3  to  the  function  Wh  and  use  inequality  (8)  to  complete 
the  proof  of  the  Theorem  1 .  We  also  note  that  there  are  no  difficulties  extending  the 
arguments  to  three  dimensions.  fl 

We  note  that  for  the  case  of  two  subregions,  it  is  easy  to  show  that  the  result  of 
this  theorem  cannot  be  improved. 
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4.  Multi-level  additive  Schwarz  methods.  The  method,  described  in  the 
two  previous  sections,  uses  two  levels  corresponding  to  a  coarse  and  a  fine  mesh  and  a 
decomposition  of  Q  into  overlapping  subregions.  If  H  is  small,  the  local  problems  are 
small  but  we  then  have  to  contend  with  a  relatively  large  linear  system  corresponding 
to  V^''  =  V" .  It  is  natural  to  try  to  decrease  the  cost  by  replacing  the  exact  solver  on 
the  coarse  mesh  by  using  the  algorithm  recursively.  This  introduces  additional  levels 
of  triangulation. 

Let  us  consider  i+l  rather  than  two  levels  of  triangulation  of  f2  with  substructures 
Q.ki  and  discretization  parameters  hk,k  =  0,-  ■  ■  ,i.  We  assume  that  the  triangulation 
on  level  k  is  a  refinement  of  that  of  level  k  —  l.  On  each  level,  we  define  a  finite  element 
space  V'"'  with  V'"  =  V''.  Each  F'"=  is  represented,  as  in  the  two- level  case,  by 

The  subspace  Vjt,  is  associated  with  a  subregion  QJ-,  obtained  by  extending  Q^i  as  in 
Section  2.  Thus,  VJt,  =  Hl(Q.'i^,)f]V'"',  extended  by  zero  outside  of  dO,',,,.  The  entire 
space  can  then  be  represented  as 

(9)  v''  =  v'">  +  J2Yl Vk, . 

k=l 1=1 

The  projections  P^,  :  V^  — >  Vi,,  are  defined  as  before  by 

and  the  original  problem  is  replaced  by 

Puh  ^{Po  +  Y.Y.  Pk^>f^  =  9h  ■ 

k=l 1=1 

The  right  hand  side  is  constructed  as  in  Section  2. 

We  note  that  we  are  now  able  to  design  methods  for  which  there  is  fixed  upper 
bound  on  the  dimension  of  the  subspaces;  when  h  decreases,  we  can  gradually  increase 
i  keeping  the  size  of  the  subproblems  uniformly  bounded.  In  this  case,  the  work 
per  iteration  step  grows  only  linearly  with  the  number  of  degrees  of  freedom  of  the 
whole  problem.  The  condition  number  of  the  stiffness  matrices  of  the  small  problems, 
corresponding  to  individual  subspaces,  are  uniformly  bounded  and  these  matrices 
can  be  replaced  by  diagonal  matrices  without  affecting  the  condition  number  of  the 
iteration  operator  by  more  than  a  constant  factor. 

Dryja  and  Widlund  [19],  have  previously  established  that  the  condition  number  of 
P  grows  at  most  quadratically  with  £,  the  number  of  levels.  Recently,  Xuejun  Zhang 
[37]  has  established  a  constant  upper  bound  for  the  eigenvalues  of  P  and  obtained: 

Theorem  4.  There  exist  constants  70  and  7^,  independent  of  the  mesh  sizes  and 
^,  the  number  of  levels,  such  that  the  operator  P  defined  by  the  decomposition  (9), 
satisfies 

7o(^  +  l)~^a(u/,,  u/i)  <  a{Puh,  Uh)  <  l\a{uh,  Uh)  ■ 


In  the  case  of  H^ —  regular Uy,   e.g.    if  the  region  is  convex,   there  is  also  a  constant 
lower  bound. 

Zhang  also  notes  that  the  case,  where  the  refinement  from  each  level  to  the  next 
is  by  a  factor  two  and  all  the  subspaces  Vk,,k  >  1,  are  of  dimension  one,  is  covered  by 
his  theory.  The  resulting  algorithm  can  also  be  derived  from  the  BPX  algorithm,  see 
Bramble,  Pasciak  and  Xu  [5],  by  using  a  simple  diagonal  scaling.  As  a  consequence, 
Zhang  has  improved  the  bound  for  the  BPX  algorithm.  Previously  only  a  CP  bound 
was  known,  in  the  absence  of  an  extra  regularity  assumption. 

5.  Nonsymmetric  and  indefinite  problems.  Let  fi  be  an  open,  bounded 
polygonal  region  in  R'^,d  =  2  or  3,  with  boundary  dO,.  Consider  the  homogeneous 
Dirichlet  boundary  value  problem: 

(10)  {  Lu    =    f    in     Q, 

^      '  \      u     =     0     on    on. 

The  elliptic  operator  L  has  the  form 

i,j=i  <^^«  oxj  fr{  dx, 

All  the  coefficients  are,  by  assumption,  sufficiently  smooth  and  the  matrix  {a,_,(T)} 
is  symmetric  and  uniformly  positive  definite  for  x  e  ^.  The  right  hand  side  /  G  L^(fi). 
We  also  assume  that  the  equation  has  a  unique  solution  in  Hq{Q,). 

The  weak  form  of  equation  (10)  is:  Find  u  G  Hq(Q)  such  that 

(11)  b{u,v)  =  f{v),     \/veH^ia). 

The  bilinear  form  b{u,v)  is  defined  by 

J^    f        du  dv  ,       J^      r      du  r 

b{u,v)   =    2^    /    a,j--—-—dx  +  }^2      b,-—vdx+  /    cuvdx 
,~^i  J  Q      ox  J  dx,  f^    Jq     dx,  J  Q 


or 


u        A  V-    /"         du   dv  ^        ^   r      du         d{ku)     ,  f 

b{u,v)   =    2^    /    a,j-—--dx  +  }^      b,—-v  +  -\ — -vdx  +  /    cuvdx. 
x%\  J  "      OX]  ox,  Jz^  Jq     dxi  dx,  J  Q 

Here,  c(a:)  =  c{x)  -  Ef=i  db,{x)/dxi. 


We  also  use  two  other  bilinear  forms 

du   dv 


I       ^         v^    /■        du  dv 
a{u,v)   =    J2        a„  ——dx 


and 

du         d{b,u) 


(       \         sr^  f  L  ^^         d(b,u]     , 
siu,v)   =   }2       h,~-v  + -^-^-^vdx 
jr[  Jq    dx,         dx. 
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which  correspond  to  the  second  order  terms  and  the  skew-symmetric  part  of  L,  re- 
spectively. The  bihnear  form  a(-,  •)  defines  a  norm,  which  we  denote  by  ||  •  ||a.  Under 
the  given  assumptions  on  the  coefficients  Oij,  this  norm  is  equivalent  to  the  Hq  norm. 
Symmetry  is  always  with  respect  to  a(-,-)  and  the  adjoint  S  of  an  operator  S  is 
defined  by  a{S^u,  v)  =  a{u,  Sv).  It  is  also  easy  to  verify  that 

s{u,v)    =    —s(v,u),    yu^v  G  Hq{Q,). 

Using  elementary,  standard  tools,  it  is  easy  to  establish  the  following  inequalities: 

(i)  I  b{u,v)  \<  C\\u\\a\\v\U,    Vu,u  G  H'oin). 

(ii)  Garding's  inequality:  There  exists  a  constant  C,  such  that 

\W\\l-C\\u\\l^^)<b(u,u),  Vueif^(fi). 

(iii)  There  exists  a  constant  C,  such  that 

I  s{u,v)  \<  C||u||a||i'||L2(n),    Vu,t;  G  Ho{9.), 

1  s(u,v)  \<  C||t;||J|u||z,2(n),    ^u.v  e  H^m. 

We  note  that  the  bounds  for  b{-,-)  and  s{-,-)  are  different,  since  each  term  in 
s{-,-)  contains  a  factor,  which  is  of  zero  order.  This  enables  us  to  control  the  skew- 
symmetric  term  and  makes  our  analysis  possible;  see  further  discussion  at  the  end  of 
this  section. 

We  also  use  the  following  regularity  result,  cf.  Grisvard  [22]  or  Necas  [27]. 

(iv)  The  solution  w  of  the  adjoint  equation 

b{6,w)   =  g{4>),    y<t>eHl{Q) 
satisfies 

||«'||//'+^(n)  <  ^11^111,2(0)  , 

where  7  depends  on  the  interior  angles  of  5fi,  is  independent  of  g  and  is  at  least  1/2. 

We  approximate  equation  (11)  by  a  Galerkin  conforming  finite  element  method 
as  in  (1).  For  simplicity,  we  consider  only  continuous,  piecewise  linear,  triangular 
elements  in  E?  and  tetrahedral  elements  in  H?. 

We  consider  the  same  additive  Schwarz  algorithm,  which  were  introduced  previ- 
ously in  this  paper,  and  also  its  multiplicative  counterpart.  Our  results  have  also  been 
established  for  several  other  classes  of  algorithms,  cf.    Cai  and  Widlund  [11].   There 
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is  also  a  more  abstract  theory,  which  we  will  describe  briefly  without  providing  any 
proofs. 

We  introduce  operators  which  can  replace  the  projections  P,  and  the  operators 
T,  used  previously.  They  are  no  longer  symmetric. 

Definition:     For  i  =  0,-  ■    ,N: 

For  any  w^  G  V^,  P,w^  G  V/"  13  the  solution  of  the  finite  element  equation 

For  any  w''  G  V",  T^w''  G  V^^  is  the  solution  of  the  finite  element  equation 

a{T,w\v^)   =  b{w\v^),    Wv^eV,'. 

Other  choices  for  the  operators  T,  are  also  possible.  We  introduce  operators  which 
define  two  possible  additive  Schwarz  methods 

P   =   Po  +  Pi   +■■■   +  Pn 
and 

T  =  Po  +  Ti  +•••  +  Tiv. 

The  only  difference  between  P  and  T  is  that,  for  i  >  0,  we  replace  the  operator 
Pi,  corresponding  to  Q[,  by  T,.  The  coarse  mesh  operator  is  not  changed;  it  should 
always  correspond  to  the  exact  solution  of  the  original  problem  in  the  coarse  space. 

The  problems  that  we  are  considering  are  always  nonsymmetric  and  we  there- 
fore have  to  replace  the  conjugate  gradient  method  by  a  more  appropriate  iterative 
method.  In  our  experiments  with  additive  methods,  we  have  always  used  the  GMRES 
method,  cf.  Saad  and  Schultz  [28],  and  Eisenstat,  Elman  and  Schultz  [21].  This  is 
a  generalized  minimum  residual  method  that  in  practice  has  proven  quite  powerful 
for  a  large  class  of  nonsymmetric  problems.  The  method  converges  satisfactorily  for 
positive  definite  problems  and  a  bound  on  the  decay  of  the  norm  of  the  residual  in 
n  iterations  is  given  by  (1    —    c^/C^Y'^  .  Here  Cp     =     inf^^o  a(a;,Tx)/a(x,a;)  and 

Cp   =   sup^_ioll^^lU/lkl|a- 

In  our  discussion  of  the  symmetric,  positive  definite  case,  we  have  found  that  the 
analysis  of  both  additive  and  multiplicative  cases  requires  estimates  of  the  constant 
Co,  the  spectral  radius  of  the  matrix  S  and  the  norms  of  the  operators  T,.  In  the 
nonsymmetric  case,  we  must  do  additional  work. 

For  the  multiphcative  algorithms,  we  wish  to  estimate  the  spectral  radius  of  the 
error  propagation  operator  Ej  given  in  formula  (7).  Following  Bramble  et  al.  [4], 
we  begin  by  observing  that,  with  Ej  =  (I  -  T,)  ■  ■  ■  [I  -  Ti)(/  -  To),  E_i  =  I,  and 
R,=T,+Tl-TJ^T,, 

Ej  Ej  -  Ej^^Ej+i  =  Ej  Rj+iEj. 
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This  leads  to  the  identity 

J 

(12)  /  -  E^Ej  =  iZo  +  E  Ej_,R,E,.r  . 

To  simpHfy  the  notations,  we  here  use  To  =  Pq. 

A  satisfactory  upper  bound  for  p(Ej)  could  therefore  be  obtained  by  showing  that 
the  operator  on  the  right  hand  side  of  (12)  is  sufficiently  positive  definite.  It  therefore 
seems  natural  to  assume  that  the  operators  Ri  axe  positive  semidefinite.  It  is  easy  to 
see  that  this  is  so  if  T^  =  T,  >  0,  but  such  an  assumption  on  R^  can  often  not  be 
established  in  our  nonsymmetric  applications.  In  the  general  case,  we  therefore  make 
a  different 

Assumption  1.  There  exist  paTameters  Si,  which  can  be  made  sufficiently  small, 
and  a  constant  7  >  0,  such  that 

(13)  R,  =  T,  +  Tj  -  TjT,  >  fT^T,  -  6,1  . 

We  note  that  if  we  can  bound  T,  +  T^  +  6,/  from  below  by  some  positive  constant 
multiple  of  T^Ti,  then  Assumption  1  is  satisfied  for  aTi  for  a  sufficiently  small  a. 
It  is  well  known  that  such  a  rescaling  often  is  necessary  to  obtain  convergence  for 
nonsymmetric  problems. 

It  follows  easily  from  Assumption  1  that 

(14)  iir.ii.  <  u,  =  (1  +  yi  +  6.(i  +  7))/(i  +  7)  <  7^  + 1 . 

1  +  72 

To  simplify  our  calculations  and  formulas,  we  always  assume  that  uj  =  max.Wi  >  1. 
It  is  easy  to  show  that  ||/  —  r,||a  <  1  +  ^  and  that  for  a  small  enough  (5,-,  ||r,||a  <  2. 
In  the  case  when  T;  is  positive  definite,  symmetric,  6i  =  0, 

(15)  0<r,  <u;/  =  2/(l+7)/ 
and 

(16)  i?.>(2-u;)T.  . 

In  our  discussion  of  the  symmetric  case  and  the  additive  algorithms,  we  have  seen 
that  we  can  often  establish  a  lower  bound  on  the  operator  T  =  To  +  •••  Tj  in  terms 
of  the  constant  Co  of  Lemma  1.  By  using  similar  tools,  we  can  often  establish 

Assumption  2.    There  exists  a  constant  Co  >  0,  such  that 

(17)  YlT^T,>Co'l. 

We  note  that  in  the  symmetric  case,  P^  -  Pi  and  generally  T,  >  T^^/2,  if  ||T,||a  <  2. 
In  this  case,  Assumption  2  is  therefore  closely  related  to  the  bound  given  in  Lemma 
1. 
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We  can  establish  Assumption  2,  for  a  sufficiently  fine  coarse  mesh,  by  using  only 
Lemma  1,  Garding's  inequality,  cf.  inequality  (ii)  above,  and  a  result  of  Schatz's 
[29].  We  note  that  Schatz  proved  the  existence  of  finite  elment  solutions  for  elliptic 
problems,  such  as  those  introduced  in  the  beginning  of  this  section,  for  a  sufficiently 
small  mesh  size.  (In  the  variant  that  we  need  here,  we  replace  the  approximate 
solution  by  the  coarse  mesh  solution  and  the  exact  solution  of  the  continuous  problem 
by  the  finite  element  solution  in  V'^.)  In  the  proof  of  the  variant  of  Schatz's  result, 
Garding's  inequality,  (ii),  and  the  regularity  result,  (iv),  but  not  (iii)  are  used.  Details 
are  given  for  the  special  Schwarz  method,  considered  in  this  paper,  in  Lemma  5  of 
Cai  and  Widlund  [11].  That  paper  also  contains  results  for  Yserentant's  method  and 
iterative  substructuring  methods  for  problems  in  the  plane. 

If  we  combine  Assumptions  1  and  2,  we  obtain  the  required  positive  lower  bound 
for  T  that  is  required  for  the  estimate  of  the  rate  of  convergence  of  the  GMRES 
algorithm  provided  that  '^6,  is  made  sufficiently  small,  cf.  Cai  and  Widlund  [10]  for 
details.  In  the  applications  this  requirement  is  met  by  decreasing  the  coarse  mesh 
size.  The  upper  bound  is  given  by 

II  E^.lla  <  2(M^)  +  1)/(1  +  7)  +  l/2i:<5.  . 

1=0  i=0 

This  bound  is  obtained  quite  straightforwardly. 

For  the  multiplicative  case  the  following  result  holds.  A  proof  will  be  given  in  Cai 
and  Widlund  [10]  together  with  several  applications. 

Theorem  5.   In  the  general  case,  the  multiplicative  algorithm  is  convergent  if 
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dominates 

J  j-i 

]=Q  :=0 

by  a  sufficiently  large  constant  factor.    Under  this  assumption,  there  exists  a  constant 
C  >  0,  such  that 


pl.Ej)  < 


i 


1  - 


c 


i^'\£\l+ {i:s,y)Q 


The  main  effort  in  the  application  of  this  theory  goes  into  the  verification  of 
Assumption  1.  So  far  we  have  have  only  been  able  to  do  so  by  using  the  inequality 
(iii)  and  a  perturbation  argument.  However,  there  are  eUiptic  problems  which  satisfy 
Garding's  inequality  but  for  which  the  skew-symmetric  part  is  not  a  relatively  compact 
perturbation.  In  such  a  case  an  inequality  serving  the  same  purpose  as  (iii)  will  not 
hold.  For  an  example  of  elUptic  problems  of  this  type,  see  Bramble,  Leyk  and  Pasciak 
[3]  in  which  several  interesting  algorithms,  which  differ  from  ours,  are  considered. 
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